2. Gene expression analysis
The output from cellranger’s count function returns gene expression information in the form of an expression matrix, barcodes, and gene identifiers. The function automate_GEX allows us to automates the transcriptional analysis of the gene expression libraries from cellranger. The input directory should be set to the directory containing these three files (barcodes.tsv.gz, features.tsv.gz, matrix.mtx.gz). If it is desired to analyze multiple distinct transcriptomes (e.g. on separate UMAP space), these directories should go into separate list objects. For example if there are 20 repertoires and it is wanted to analyze them separately (20 separate UMAPs, 20 separate repertoires) then each outer list element of the input directories will contain the directory for the individual repertoire/GEX files. The output of this function is a Seurat object similar to the standard pipeline demonstrated in the Seurat vignette. This involves scaling, normalizing, clustering, and performing dimension reductionality (tSNE and UMAP by default). At this stage there is no incorporation of repertoire features, which need to be explicitly integrated using subsequent functions in Platypus.
In the example below two directories are supplied to two distinct list elements, and thus these distinct transcriptomes (corresponding to two distinct repertoires) will each be associated to a distinct Seurat object.
### Removing any previous versions of the package
# First we will ensure that there is no previous version installed locally
#detach("package:Platypus", unload=TRUE)
#remove.packages("Platypus")
### Downloading and installing Platypus
# First we need to download the most recent version from the master branch at https://github.com/alexyermanos/Platypus we can install the package using the following command.
# WARNING: This needs to be replaced with your own directory where the downloaded package is found
# For MacOS useres it may look like this
#install.packages("~/Downloads/Platypus_2.0.4.tar.gz", repos = NULL, type="source")
# For windows it will likely look something like this.
# WARNING: You will need to replace YourPCName with your User name for the windows account in the directory.
# install.packages("C:\Users\YourPCName\Downloads\Platypus_2.0.4.tar.gz", repos = NULL, type="source")
# Now we can load the installed package into the R environment. In case of problems with instally other R packages that are used in Platypus, please see the README file at the https://github.com/alexyermanos/Platypus, where we outline how to install the other R packages for both Windows and MacOS.
#library(Platypus)
# The individual R functions can additionally be found on the github in the Functions branch. Within this branch, there is an folder "R" which contains the individual functions. This can similarly be downloaded and loaded into the R environment incase not all functions are desired. Similarly, these functions are actively updated and may include more features than the in original tar.gz file.
### Downloading the test data
# The COVID-19 data (~136 MB size of the zip file) can be found at the following link https://polybox.ethz.ch/index.php/s/fxQJ3NrRSwiPSSo This dataset contains VDJ (separate libraries for B and T cells) and GEX libraries from two convalescent COVID-19 patients.
# After downloading the zip file named "PlatypusTestData.zip", please unzip the file and find the path to the newly formed folder. Typically this will be in the Downloads folder, so the code below should work on MacOS. For windows please uncomment the code and change the user name to match your PC.
directory_to_covid_patients_gex <- list()
directory_to_covid_patients_gex[[1]] <- c("~/Downloads/PlatypusTestData/Patient1_GEX/")
directory_to_covid_patients_gex[[2]] <- c("~/Downloads/PlatypusTestData/Patient2_GEX/")
# For windows:
#directory_to_covid_patients_gex[[1]] <- c("C:\Users\YourPCName\Downloads\PlatypusTestData\Patient1_GEX")
#directory_to_covid_patients_gex[[2]] <- c("C:\Users\YourPCName\Downloads\PlatypusTestData\Patient2_GEX")GEX_automate allows us to perform the standard Seurat pipeline in a single line of code with the ability to change the parameters used by Seurat involving minimum read numbers, mitochondrial gene percentages, cluster resolution, etc. In addition we can remove antibody and TCR genes from the data set (e.g. IGHV1-1), to not let clonality impact transcriptional clustering, by simply setting the VDJ.gene.filter argument to TRUE.
## [1] "This may take longer than other repertoire associated functions. Please see Seurat vingenettes for further information"
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 6871
## Number of edges: 227893
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8914
## Number of communities: 13
## Elapsed time: 0 seconds## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 6055
## Number of edges: 205942
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8829
## Number of communities: 12
## Elapsed time: 0 secondsBased on the input given (2 directories as separate list objects), we get 2 distinct Seurat objects.
## [1] 2## [1] "list"## [1] "Seurat"
## attr(,"package")
## [1] "Seurat"## [1] 6871## [1] 6055Now we can visualize the 2D plots for each patient individually. By default, UMAP, PCA, and TSNE reductions are included in the object. Under default parameters this will display all cells in the GEX library from Patient1. By changing to the second list element we can view the cells from Patient2.
Seurat::DimPlot(covid_gex_patients_not_integrated[[1]],reduction = "umap")
Seurat::DimPlot(covid_gex_patients_not_integrated[[1]],reduction = "pca")
Seurat::DimPlot(covid_gex_patients_not_integrated[[1]],reduction = "tsne")
# UMAP for the second patient
Seurat::DimPlot(covid_gex_patients_not_integrated[[2]],reduction = "umap")
If we are intrested in merging the transcripotmes (corresponding to two distinct gene expression libraries), then the two directories are supplied to a single list element. The output from GEX_automate (corresponding to a single Seurat object list element this time) contains a column with the sample origin infromation based on the order of the input directories. Additionally, one can label different samples by group in the case that biological replicates are present. In this case we simply set patients to group1 and group2
directory_to_covid_patients_integrated <- list()
directory_to_covid_patients_integrated[[1]] <- c("~/Downloads/PlatypusTestData/Patient1_GEX/",
"~/Downloads/PlatypusTestData/Patient2_GEX/")
## Here we use the previous version of automate_GEX to produce the identical UMAP seen in the manuscript.
covid_gex <- Platypus::automate_GEX(GEX.outs.directory.list = directory_to_covid_patients_integrated[1:1],integration.method = "scale.data",mito.filter = 20,cluster.resolution = 0.5,VDJ.gene.filter = T)## [1] "This may take longer than other repertoire associated functions. Please see Seurat vignettes for further information"
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 12926
## Number of edges: 426416
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8933
## Number of communities: 12
## Elapsed time: 2 seconds### In this case, patient1 is the covid_gex[[1]]$sample_id==1 and patient2 is the covid_gex[[1]]$sample_id==2.
class(covid_gex) ## list## [1] "list"## [1] 1## [1] "Seurat"
## attr(,"package")
## [1] "Seurat"We can extract which cells come from which sample based on the sample_id in the Seurat object.
##
## 1 2
## 6871 6055We can see that the first patient had 6871 cells and the second patient had 6055 cells found in their gene expression data based on the current filtering.
Now we can visualize the samples using a UMAP/tsne plot.
We observe that under this Seurat pipeline there are 12 distinct clusters. The number of clusters can be changed by altering the cluster.resolution argument in the automate_GEX function.
In this current UMAP we do not see which cells are coming from which patient. To see this, we can simply separate the cells in the DimPlot by the sample_id vector in the Seurat object.
Seurat::DimPlot(covid_gex[[1]],reduction = "umap",split.by = "sample_id")
Seurat::DimPlot(covid_gex[[1]],reduction = "umap")
We can visually observe that the majority of clusters have cells from both patients, suggesting a similar distribution of transcriptional properties between the two samples. We can also plot the cluster membership for each of the distinct samples (patients) using the GEX_cluster_per_sample function.
Finally we can also look at the B, CD4 and CD8 clusters for each patient.
We can see that both patients seem to have B cells (CD19+), CD4 T cells, and CD8 T cells. Platypus also allows us to assigning cell type/state identity to the different clusters by using the GEX_phenotype. This functon takes the Seurat object as input and uses canonical markers to easily match the clustering to known cell types. The user also has the possibility to use a custom list of markers and their associate cell types/states.
## Loading required package: do## Warning: package 'do' was built under R version 3.6.2##
## Attaching package: 'do'## The following object is masked from 'package:Biostrings':
##
## reverse## The following object is masked from 'package:XVector':
##
## reverse## The following objects are masked from 'package:IRanges':
##
## mid, reverse## The following object is masked from 'package:S4Vectors':
##
## expand## Loading required package: useful##
## Attaching package: 'useful'## The following objects are masked from 'package:do':
##
## left, rightcovid_gex[[1]] <- Platypus::GEX_phenotype(covid_gex[[1]], default = F,
cell.state.markers=c("CD8A+;CCL5+;CD44+;IL7R-;CD19-",
"CD8A+;CCL5-;CD44+;IL7R+;CD19-"),
cell.state.names=c("EffectorCD8",
"MemoryCD8"))The resulting Seurat object now contains a “cell.state” column which can be used for annotation in the DimPlot function of the Seurat package.
